QTS data: Spiritual curiosity relationships

Last updated 2018-04-12

This is a first-pass exploration of some of the relationships among QTS data for the “Spiritual curiosity” interview.

Time praying & spiritual experiences

Requested by Josh

Time praying as indexed by prayfreqmin

How many minutes each day do you do that on average throughout the week? [use real world examples, ask alone? With others? Probe away– Answer with interviewer judgement for answers]

prayfreqmin_count <- d %>%
  filter(!is.na(prayfreqmin), !is.na(godvoxaloud),
         prayfreqmin != "Other", godvoxaloud != "Other") %>%
  count(country, prayfreqmin) %>%
  data.frame()
prayfreqmin_count_by_quad <- d %>%
  filter(!is.na(prayfreqmin), !is.na(godvoxaloud),
         prayfreqmin != "Other", godvoxaloud != "Other") %>%
  count(country, urban_rural, charismatic_local, prayfreqmin) %>%
  data.frame()

Spiritual experiences indexed by godvoxaloud

Some people say that they have heard God* speak out loud to them. Has this ever happened to you?

d %>%
  filter(!is.na(prayfreqmin), !is.na(godvoxaloud)) %>% #,
         # prayfreqmin != "Other", godvoxaloud != "Other") %>%
  ggplot(aes(x = prayfreqmin, y = godvoxaloud, color = researcher)) +
  facet_grid(urban_rural ~ charismatic_local ~ country) +
  geom_point(alpha = 0.3, position = position_jitter(0.2, 0.2)) +
  scale_color_manual(values = custom_pal) +
  labs(title = paste("prayfreqmin", "godvoxaloud", sep = " x "),
       # subtitle = "Excluding people who did not have a clear answer",
       x = "prayfreqmin",
       y = "godvoxaloud",
       color = "Researcher") +
  theme_bw() +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  guides(color = guide_legend(ncol = 6, byrow = TRUE, 
                              override.aes = list(alpha = 1)))

d %>%
  filter(!is.na(prayfreqmin), !is.na(godvoxaloud),
         prayfreqmin != "Other", godvoxaloud != "Other") %>%
  ggplot(aes(x = prayfreqmin, alpha = godvoxaloud, fill = researcher)) +
  facet_grid(urban_rural ~ charismatic_local ~ country) +
  geom_bar(position = "fill") +
  geom_text(data = prayfreqmin_count_by_quad,
            aes(x = prayfreqmin, y = 1, alpha = NULL, fill = NULL,
                label = paste0("(n=", n, ")")),
            size = 2, nudge_y = 0.05) +
  # scale_fill_brewer(guide = NULL, palette = "Dark2") +
  scale_fill_manual(guide = NULL, values = custom_pal) +
  scale_alpha_discrete(range = c(0.5, 1)) +
  scale_y_continuous(breaks = seq(0, 1, 0.25)) +
  labs(title = paste("prayfreqmin", "godvoxaloud", sep = " x "),
       subtitle = "Excluding people who did not have a clear answer",
       x = "prayfreqmin",
       y = "Proportion",
       alpha = "godvoxaloud",
       fill = "Researcher") +
  theme_bw() +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  guides(alpha = guide_legend(),
         fill = "none")

Ignoring unknown aesthetics: fill

d %>%
  filter(!is.na(prayfreqmin), !is.na(godvoxaloud),
         prayfreqmin != "Other", godvoxaloud != "Other") %>%
  ggplot(aes(x = prayfreqmin, alpha = godvoxaloud, fill = country)) +
  facet_grid(. ~ country) +
  geom_bar(position = "fill") +
  geom_text(data = prayfreqmin_count,
            aes(x = prayfreqmin, y = 1, alpha = NULL, fill = NULL,
                label = paste0("(n=", n, ")")),
            size = 2, nudge_y = 0.05) +
  scale_fill_brewer(guide = NULL, palette = "Dark2") +
  # scale_fill_manual(guide = NULL, values = custom_pal) +
  scale_alpha_discrete(range = c(0.5, 1)) +
  scale_y_continuous(breaks = seq(0, 1, 0.25)) +
  labs(title = paste("prayfreqmin", "godvoxaloud", sep = " x "),
       subtitle = "Excluding people who did not have a clear answer",
       x = "prayfreqmin",
       y = "Proportion",
       alpha = "godvoxaloud",
       fill = "Researcher") +
  theme_bw() +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  guides(alpha = guide_legend(),
         fill = "none")

Ignoring unknown aesthetics: fill

d1 <- d %>%
  filter(!is.na(prayfreqmin), !is.na(godvoxaloud),
         prayfreqmin != "Other", godvoxaloud != "Other") %>%
  mutate(prayfreqmin_num = as.numeric(prayfreqmin) - 1,
         godvoxaloud = factor(godvoxaloud,
                               levels = c("No", "Yes")),
         godvoxaloud_num = as.numeric(godvoxaloud) - 1) %>%
  select(country, researcher, urban_rural, charismatic_local, subject_name,
         starts_with("prayfreqmin"), starts_with("godvoxaloud"))
contrasts(d1$country) <- cbind("GH" = c(-1, 1, 0, 0, 0),
                               "TH" = c(-1, 0, 1, 0, 0),
                               "CH" = c(-1, 0, 0, 1, 0),
                               "VT" = c(-1, 0, 0, 0, 1))
contrasts(d1$prayfreqmin) <- contr.poly(4)
r1 <- glm(godvoxaloud_num ~ prayfreqmin + country,
                  family = "binomial", data = d1)
summary(r1)

Call:
glm(formula = godvoxaloud_num ~ prayfreqmin + country, family = "binomial", 
    data = d1)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.703  -1.027  -0.455   1.032   2.154  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)   -0.47926    0.13784  -3.477 0.000507 ***
prayfreqmin.L  0.73747    0.24453   3.016 0.002562 ** 
prayfreqmin.Q  0.54372    0.27178   2.001 0.045436 *  
prayfreqmin.C -0.35794    0.29205  -1.226 0.220347    
countryGH      0.97536    0.21961   4.441 8.94e-06 ***
countryTH      0.25817    0.26628   0.970 0.332274    
countryCH      0.30449    0.24600   1.238 0.215809    
countryVT      0.05578    0.25920   0.215 0.829609    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 446.17  on 328  degrees of freedom
Residual deviance: 385.93  on 321  degrees of freedom
AIC: 401.93

Number of Fisher Scoring iterations: 4

Spiritual experiences indexed by godviavisions

_Some people say that they have had a vision from God*—they have a picture, but it is like they see it with their eyes. Has anything like that happened to you?_

d %>%
  filter(!is.na(prayfreqmin), !is.na(godviavisions)) %>% #,
         # prayfreqmin != "Other", godviavisions != "Other") %>%
  ggplot(aes(x = prayfreqmin, y = godviavisions, color = researcher)) +
  facet_grid(urban_rural ~ charismatic_local ~ country) +
  geom_point(alpha = 0.3, position = position_jitter(0.2, 0.2)) +
  scale_color_manual(values = custom_pal) +
  labs(title = paste("prayfreqmin", "godviavisions", sep = " x "),
       # subtitle = "Excluding people who did not have a clear answer",
       x = "prayfreqmin",
       y = "godviavisions",
       color = "Researcher") +
  theme_bw() +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  guides(color = guide_legend(ncol = 6, byrow = TRUE, 
                              override.aes = list(alpha = 1)))

d %>%
  filter(!is.na(prayfreqmin), !is.na(godviavisions),
         prayfreqmin != "Other", godviavisions != "Other") %>%
  ggplot(aes(x = prayfreqmin, alpha = godviavisions, fill = researcher)) +
  facet_grid(urban_rural ~ charismatic_local ~ country) +
  geom_bar(position = "fill") +
  geom_text(data = prayfreqmin_count_by_quad,
            aes(x = prayfreqmin, y = 1, alpha = NULL, fill = NULL,
                label = paste0("(n=", n, ")")),
            size = 2, nudge_y = 0.05) +
  # scale_fill_brewer(guide = NULL, palette = "Dark2") +
  scale_fill_manual(guide = NULL, values = custom_pal) +
  scale_alpha_discrete(range = c(0.5, 1)) +
  scale_y_continuous(breaks = seq(0, 1, 0.25)) +
  labs(title = paste("prayfreqmin", "godviavisions", sep = " x "),
       subtitle = "Excluding people who did not have a clear answer",
       x = "prayfreqmin",
       y = "Proportion",
       alpha = "godviavisions",
       fill = "Researcher") +
  theme_bw() +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  guides(alpha = guide_legend(),
         fill = "none")

Ignoring unknown aesthetics: fill

d %>%
  filter(!is.na(prayfreqmin), !is.na(godviavisions),
         prayfreqmin != "Other", godviavisions != "Other") %>%
  ggplot(aes(x = prayfreqmin, alpha = godviavisions, fill = country)) +
  facet_grid(. ~ country) +
  geom_bar(position = "fill") +
  geom_text(data = prayfreqmin_count,
            aes(x = prayfreqmin, y = 1, alpha = NULL, fill = NULL,
                label = paste0("(n=", n, ")")),
            size = 2, nudge_y = 0.05) +
  scale_fill_brewer(guide = NULL, palette = "Dark2") +
  # scale_fill_manual(guide = NULL, values = custom_pal) +
  scale_alpha_discrete(range = c(0.5, 1)) +
  scale_y_continuous(breaks = seq(0, 1, 0.25)) +
  labs(title = paste("prayfreqmin", "godviavisions", sep = " x "),
       subtitle = "Excluding people who did not have a clear answer",
       x = "prayfreqmin",
       y = "Proportion",
       alpha = "godviavisions",
       fill = "Researcher") +
  theme_bw() +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  guides(alpha = guide_legend(),
         fill = "none")

Ignoring unknown aesthetics: fill

d2 <- d %>%
  filter(!is.na(prayfreqmin), !is.na(godviavisions),
         prayfreqmin != "Other", godviavisions != "Other") %>%
  mutate(prayfreqmin_num = as.numeric(prayfreqmin) - 1,
         godviavisions = factor(godviavisions,
                               levels = c("No", "Yes")),
         godviavisions_num = as.numeric(godviavisions) - 1) %>%
  select(country, researcher, urban_rural, charismatic_local, subject_name,
         starts_with("prayfreqmin"), starts_with("godviavisions"))
contrasts(d2$country) <- cbind("GH" = c(-1, 1, 0, 0, 0),
                               "TH" = c(-1, 0, 1, 0, 0),
                               "CH" = c(-1, 0, 0, 1, 0),
                               "VT" = c(-1, 0, 0, 0, 1))
contrasts(d2$prayfreqmin) <- contr.poly(4)
r2 <- glm(godviavisions_num ~ prayfreqmin + country,
                  family = "binomial", data = d2)
summary(r2)

Call:
glm(formula = godviavisions_num ~ prayfreqmin + country, family = "binomial", 
    data = d2)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.8640  -1.0136  -0.5422   1.1679   2.0691  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)    -0.6210     0.1361  -4.562 5.06e-06 ***
prayfreqmin.L   0.9386     0.2470   3.800 0.000145 ***
prayfreqmin.Q   0.8468     0.2739   3.092 0.001989 ** 
prayfreqmin.C   0.3692     0.2926   1.262 0.207031    
countryGH       1.0290     0.2239   4.597 4.29e-06 ***
countryTH       0.6839     0.2719   2.515 0.011894 *  
countryCH      -0.9077     0.2933  -3.095 0.001968 ** 
countryVT       0.1282     0.2570   0.499 0.617930    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 443.63  on 334  degrees of freedom
Residual deviance: 386.58  on 327  degrees of freedom
AIC: 402.58

Number of Fisher Scoring iterations: 4

Spiritual experiences indexed by godviabodyexperiences

Some people have particular experiences in your body that they associate with God* or spirit. Does that happen for you? [examples: warm hands, goosebumps, fire in the belly]

d %>%
  filter(!is.na(prayfreqmin), !is.na(godviabodyexperiences)) %>% #,
         # prayfreqmin != "Other", godviabodyexperiences != "Other") %>%
  ggplot(aes(x = prayfreqmin, y = godviabodyexperiences, color = researcher)) +
  facet_grid(urban_rural ~ charismatic_local ~ country) +
  geom_point(alpha = 0.3, position = position_jitter(0.2, 0.2)) +
  scale_color_manual(values = custom_pal) +
  labs(title = paste("prayfreqmin", "godviabodyexperiences", sep = " x "),
       # subtitle = "Excluding people who did not have a clear answer",
       x = "prayfreqmin",
       y = "godviabodyexperiences",
       color = "Researcher") +
  theme_bw() +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  guides(color = guide_legend(ncol = 6, byrow = TRUE, 
                              override.aes = list(alpha = 1)))

d %>%
  filter(!is.na(prayfreqmin), !is.na(godviabodyexperiences),
         prayfreqmin != "Other", godviabodyexperiences != "Other") %>%
  ggplot(aes(x = prayfreqmin, alpha = godviabodyexperiences, fill = researcher)) +
  facet_grid(urban_rural ~ charismatic_local ~ country) +
  geom_bar(position = "fill") +
  geom_text(data = prayfreqmin_count_by_quad,
            aes(x = prayfreqmin, y = 1, alpha = NULL, fill = NULL,
                label = paste0("(n=", n, ")")),
            size = 2, nudge_y = 0.05) +
  # scale_fill_brewer(guide = NULL, palette = "Dark2") +
  scale_fill_manual(guide = NULL, values = custom_pal) +
  scale_alpha_discrete(range = c(0.5, 1)) +
  scale_y_continuous(breaks = seq(0, 1, 0.25)) +
  labs(title = paste("prayfreqmin", "godviabodyexperiences", sep = " x "),
       subtitle = "Excluding people who did not have a clear answer",
       x = "prayfreqmin",
       y = "Proportion",
       alpha = "godviabodyexperiences",
       fill = "Researcher") +
  theme_bw() +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  guides(alpha = guide_legend(),
         fill = "none")

Ignoring unknown aesthetics: fill

d %>%
  filter(!is.na(prayfreqmin), !is.na(godviabodyexperiences),
         prayfreqmin != "Other", godviabodyexperiences != "Other") %>%
  ggplot(aes(x = prayfreqmin, alpha = godviabodyexperiences, fill = country)) +
  facet_grid(. ~ country) +
  geom_bar(position = "fill") +
  geom_text(data = prayfreqmin_count,
            aes(x = prayfreqmin, y = 1, alpha = NULL, fill = NULL,
                label = paste0("(n=", n, ")")),
            size = 2, nudge_y = 0.05) +
  scale_fill_brewer(guide = NULL, palette = "Dark2") +
  # scale_fill_manual(guide = NULL, values = custom_pal) +
  scale_alpha_discrete(range = c(0.5, 1)) +
  scale_y_continuous(breaks = seq(0, 1, 0.25)) +
  labs(title = paste("prayfreqmin", "godviabodyexperiences", sep = " x "),
       subtitle = "Excluding people who did not have a clear answer",
       x = "prayfreqmin",
       y = "Proportion",
       alpha = "godviabodyexperiences",
       fill = "Researcher") +
  theme_bw() +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  guides(alpha = guide_legend(),
         fill = "none")

Ignoring unknown aesthetics: fill

d3 <- d %>%
  filter(!is.na(prayfreqmin), !is.na(godviabodyexperiences),
         prayfreqmin != "Other", godviabodyexperiences != "Other") %>%
  mutate(prayfreqmin_num = as.numeric(prayfreqmin) - 1,
         godviabodyexperiences = factor(godviabodyexperiences,
                               levels = c("No", "Yes")),
         godviabodyexperiences_num = as.numeric(godviabodyexperiences) - 1) %>%
  select(country, researcher, urban_rural, charismatic_local, subject_name,
         starts_with("prayfreqmin"), starts_with("godviabodyexperiences"))
contrasts(d3$country) <- cbind("GH" = c(-1, 1, 0, 0, 0),
                               "TH" = c(-1, 0, 1, 0, 0),
                               "CH" = c(-1, 0, 0, 1, 0),
                               "VT" = c(-1, 0, 0, 0, 1))
contrasts(d3$prayfreqmin) <- contr.poly(4)
r3 <- glm(godviabodyexperiences_num ~ prayfreqmin + country,
                  family = "binomial", data = d3)
summary(r3)

Call:
glm(formula = godviabodyexperiences_num ~ prayfreqmin + country, 
    family = "binomial", data = d3)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.0197  -1.1520   0.6754   0.8918   1.3475  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)     0.7061     0.1332   5.302 1.14e-07 ***
prayfreqmin.L   0.5249     0.2378   2.207  0.02731 *  
prayfreqmin.Q   0.5157     0.2622   1.966  0.04926 *  
prayfreqmin.C   0.2188     0.2779   0.788  0.43096    
countryGH       0.5353     0.2375   2.254  0.02418 *  
countryTH       0.8843     0.3073   2.877  0.00401 ** 
countryCH      -0.4790     0.2488  -1.926  0.05415 .  
countryVT      -0.1303     0.2500  -0.521  0.60225    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 432.22  on 334  degrees of freedom
Residual deviance: 402.36  on 327  degrees of freedom
AIC: 418.36

Number of Fisher Scoring iterations: 4

Aggregate spiritual experiences score

d_spirit %>%
  filter(!is.na(prayfreqmin), !is.na(spex_score)) %>%
  distinct(country, researcher, urban_rural, charismatic_local,
           subject_name, prayfreqmin, spex_score) %>%
  ggplot(aes(x = prayfreqmin, y = spex_score, color = researcher)) +
  facet_grid(urban_rural ~ charismatic_local ~ country) +
  geom_point(alpha = 0.3, position = position_jitter(0.2, 0)) +
  scale_color_manual(values = custom_pal) +
  labs(title = paste("prayfreqmin", "spex_score", sep = " x "),
       # subtitle = "Excluding people who did not have a clear answer",
       x = "prayfreqmin",
       y = "spex_score",
       color = "Researcher") +
  theme_bw() +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  guides(color = guide_legend(ncol = 6, byrow = TRUE, 
                              override.aes = list(alpha = 1)))

d_spirit %>%
  filter(!is.na(prayfreqmin), !is.na(spex_score)) %>%
  distinct(country, researcher, urban_rural, charismatic_local,
           subject_name, prayfreqmin, spex_score) %>%
  ggplot(aes(x = prayfreqmin, y = spex_score, 
             color = country, fill = country, group = country)) +
  facet_grid(. ~ country) +
  geom_point(alpha = 0.3, position = position_jitter(0.2, 0)) +
  geom_smooth(aes(x = as.numeric(prayfreqmin)), method = "lm") +
  # geom_smooth(aes(x = as.numeric(prayfreqmin)), method = "loess") +
  # scale_color_manual(values = custom_pal) +
  # scale_fill_manual(values = custom_pal) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  labs(title = paste("prayfreqmin", "spex_score", sep = " x "),
       # subtitle = "Excluding people who did not have a clear answer",
       x = "prayfreqmin",
       y = "spex_score",
       color = "Researcher", fill = "Researcher") +
  theme_bw() +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  guides(color = guide_legend(ncol = 6, byrow = TRUE, 
                              override.aes = list(alpha = 1)),
         fill = guide_legend(ncol = 6, byrow = TRUE))

d_agg1 <- d_spirit %>%
  filter(!is.na(prayfreqmin), !is.na(spex_score),
         prayfreqmin != "Other", spex_score != "Other") %>%
  # distinct(country, researcher, urban_rural, charismatic_local,
  #          subject_name, prayfreqmin, spex_score) %>%
  mutate(prayfreqmin_num = as.numeric(prayfreqmin) - 1) # %>%
  # select(country, researcher, urban_rural, charismatic_local, subject_name,
  #        starts_with("prayfreqmin"), starts_with("spex_score"))
contrasts(d_agg1$country) <- cbind("GH" = c(-1, 1, 0, 0, 0),
                               "TH" = c(-1, 0, 1, 0, 0),
                               "CH" = c(-1, 0, 0, 1, 0),
                               "VT" = c(-1, 0, 0, 0, 1))
contrasts(d_agg1$prayfreqmin) <- contr.poly(4)
r_agg1 <- lmer(response ~ prayfreqmin + country 
               + (1 | subject_name) + (1 | question), 
               data = d_agg1)
summary(r_agg1)

Linear mixed model fit by REML ['lmerMod']
Formula: response ~ prayfreqmin + country + (1 | subject_name) + (1 |  
    question)
   Data: d_agg1

REML criterion at convergence: 9516

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.6638 -0.7863 -0.1873  0.8822  2.5579 

Random effects:
 Groups       Name        Variance Std.Dev.
 subject_name (Intercept) 0.02090  0.1446  
 question     (Intercept) 0.03517  0.1875  
 Residual                 0.17972  0.4239  
Number of obs: 7961, groups:  subject_name, 344; question, 24

Fixed effects:
               Estimate Std. Error t value
(Intercept)    0.459403   0.039550  11.616
prayfreqmin.L  0.117263   0.017788   6.592
prayfreqmin.Q  0.033800   0.019728   1.713
prayfreqmin.C  0.048443   0.021485   2.255
countryGH      0.024153   0.016975   1.423
countryTH      0.001581   0.020844   0.076
countryCH     -0.086277   0.019517  -4.421
countryVT      0.155426   0.019417   8.005

Correlation of Fixed Effects:
            (Intr) pryf.L pryf.Q pryf.C cntrGH cntrTH cntrCH
prayfrqmn.L  0.035                                          
prayfrqmn.Q -0.066  0.174                                   
prayfrqmn.C  0.024 -0.164 -0.026                            
countryGH   -0.045 -0.099  0.114  0.118                     
countryTH    0.034  0.240  0.040 -0.020 -0.259              
countryCH    0.014 -0.076 -0.095 -0.034 -0.227 -0.305       
countryVT    0.004 -0.045  0.058 -0.022 -0.209 -0.299 -0.278

r_agg2 <- lmer(response ~ prayfreqmin * country 
               + (1 | subject_name) + (1 | question), 
               data = d_agg1)
summary(r_agg2)

Linear mixed model fit by REML ['lmerMod']
Formula: response ~ prayfreqmin * country + (1 | subject_name) + (1 |  
    question)
   Data: d_agg1

REML criterion at convergence: 9563.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.6672 -0.7887 -0.1796  0.8787  2.5708 

Random effects:
 Groups       Name        Variance Std.Dev.
 subject_name (Intercept) 0.02118  0.1455  
 question     (Intercept) 0.03517  0.1875  
 Residual                 0.17973  0.4239  
Number of obs: 7961, groups:  subject_name, 344; question, 24

Fixed effects:
                         Estimate Std. Error t value
(Intercept)              0.458378   0.039986  11.463
prayfreqmin.L            0.108947   0.022863   4.765
prayfreqmin.Q            0.026234   0.023102   1.136
prayfreqmin.C            0.042343   0.023338   1.814
countryGH                0.027652   0.018607   1.486
countryTH               -0.012737   0.030968  -0.411
countryCH               -0.075852   0.022478  -3.374
countryVT                0.153817   0.020838   7.381
prayfreqmin.L:countryGH  0.006461   0.035854   0.180
prayfreqmin.Q:countryGH  0.003627   0.037213   0.097
prayfreqmin.C:countryGH  0.023226   0.038524   0.603
prayfreqmin.L:countryTH -0.039139   0.069967  -0.559
prayfreqmin.Q:countryTH -0.016074   0.061933  -0.260
prayfreqmin.C:countryTH  0.009706   0.052689   0.184
prayfreqmin.L:countryCH -0.029814   0.038224  -0.780
prayfreqmin.Q:countryCH -0.039246   0.044956  -0.873
prayfreqmin.C:countryCH  0.018775   0.050805   0.370
prayfreqmin.L:countryVT -0.012807   0.039717  -0.322
prayfreqmin.Q:countryVT  0.036971   0.041676   0.887
prayfreqmin.C:countryVT  0.035715   0.043548   0.820

Correlation matrix not shown by default, as p = 20 > 12.
Use print(x, correlation=TRUE)  or
     vcov(x)     if you need it

anova(r_agg1, r_agg2)

refitting model(s) with ML (instead of REML)

Data: d_agg1
Models:
r_agg1: response ~ prayfreqmin + country + (1 | subject_name) + (1 | 
r_agg1:     question)
r_agg2: response ~ prayfreqmin * country + (1 | subject_name) + (1 | 
r_agg2:     question)
       Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
r_agg1 11 9489.8 9566.6 -4733.9   9467.8                        
r_agg2 23 9504.7 9665.3 -4729.3   9458.7 9.134     12     0.6914

d_r_agg1_predicted <- d_agg1 %>%
  mutate(response_pred = predict(r_agg1, d_agg1)) %>%
  group_by(subject_name) %>%
  mutate(spex_score_pred = sum(response_pred, na.rm = T)) %>%
  ungroup() %>%
  group_by(country, prayfreqmin) %>%
  do(data.frame(rbind(smean.cl.boot(.$spex_score_pred))))
d_spirit %>%
  filter(!is.na(prayfreqmin), !is.na(spex_score)) %>%
  distinct(country, researcher, urban_rural, charismatic_local,
           subject_name, prayfreqmin, spex_score) %>%
  ggplot(aes(x = prayfreqmin, y = spex_score, 
             color = country, fill = country, group = country)) +
  facet_grid(. ~ country) +
  geom_point(alpha = 0.3, position = position_jitter(0.2, 0)) +
  geom_line(data = d_r_agg1_predicted, aes(y = Mean)) +
  geom_ribbon(data = d_r_agg1_predicted, 
              aes(ymin = Lower, ymax = Upper, y = NULL), 
              alpha = 0.5, size = 0) +
  # scale_color_manual(values = custom_pal) +
  # scale_fill_manual(values = custom_pal) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  labs(title = paste("prayfreqmin", "spex_score", sep = " x "),
       subtitle = "Lines are predictions from a mixed effects linear regression:\nlmer(response ~ prayfreqmin + country + (1 | subject_name) + (1 | question)",
       x = "prayfreqmin",
       y = "spex_score",
       color = "Site", fill = "Site") +
  theme_bw() +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  guides(color = "none",
         fill = "none")

Ignoring unknown aesthetics: y

Looking just at ‘apples’ (urban charismatics)

d_spirit %>%
  filter(!is.na(prayfreqmin), !is.na(spex_score),
         quad == "urban charismatic") %>%
  distinct(country, researcher, charismatic_local, urban_rural, 
           subject_name, prayfreqmin, spex_score) %>%
  ggplot(aes(x = prayfreqmin, y = spex_score, 
             color = country, fill = country, group = country)) +
  facet_grid(charismatic_local ~ urban_rural ~ country) +
  geom_point(alpha = 0.3, position = position_jitter(0.2, 0)) +
  geom_smooth(method = "lm") +
  # scale_color_manual(values = custom_pal) +
  # scale_fill_manual(values = custom_pal) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  labs(title = paste("prayfreqmin", "spex_score", sep = " x "),
       # subtitle = "Lines are predictions from a mixed effects linear regression:\nlmer(response ~ prayfreqmin + country + (1 | subject_name) + (1 | question)",
       x = "prayfreqmin",
       y = "spex_score",
       color = "Site", fill = "Site") +
  theme_bw() +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  guides(color = "none",
         fill = "none")

Looking at all charismatics

d_spirit %>%
  filter(!is.na(prayfreqmin), !is.na(spex_score),
         charismatic_local == "charismatic") %>%
  distinct(country, researcher, charismatic_local, urban_rural, 
           subject_name, prayfreqmin, spex_score) %>%
  ggplot(aes(x = prayfreqmin, y = spex_score, 
             color = country, fill = country, group = country)) +
  facet_grid(charismatic_local ~ country) +
  geom_point(alpha = 0.3, position = position_jitter(0.2, 0)) +
  geom_smooth(method = "lm") +
  # scale_color_manual(values = custom_pal) +
  # scale_fill_manual(values = custom_pal) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  labs(title = paste("prayfreqmin", "spex_score", sep = " x "),
       # subtitle = "Lines are predictions from a mixed effects linear regression:\nlmer(response ~ prayfreqmin + country + (1 | subject_name) + (1 | question)",
       x = "prayfreqmin",
       y = "spex_score",
       color = "Site", fill = "Site") +
  theme_bw() +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  guides(color = "none",
         fill = "none")

Doubt & external sensory experiences

Requested by Emily, John

“Doubt” indxed by selfunsuregodreal

Has there been a time when you yourself wondered whether God* was real?

selfunsuregodreal_count <- d %>%
  filter(!is.na(selfunsuregodreal), !is.na(godvoxaloud),
         selfunsuregodreal != "Other", godvoxaloud != "Other") %>%
  count(country, selfunsuregodreal) %>%
  data.frame()
selfunsuregodreal_count_by_quad <- d %>%
  filter(!is.na(selfunsuregodreal), !is.na(godvoxaloud),
         selfunsuregodreal != "Other", godvoxaloud != "Other") %>%
  count(country, urban_rural, charismatic_local, selfunsuregodreal) %>%
  data.frame()

Spiritual experiences indexed by godvoxaloud

Some people say that they have heard God* speak out loud to them. Has this ever happened to you?

d %>%
  filter(!is.na(selfunsuregodreal), !is.na(godvoxaloud)) %>% #,
         # selfunsuregodreal != "Other", godvoxaloud != "Other") %>%
  ggplot(aes(x = selfunsuregodreal, y = godvoxaloud, color = researcher)) +
  facet_grid(urban_rural ~ charismatic_local ~ country) +
  geom_point(alpha = 0.3, position = position_jitter(0.2, 0.2)) +
  scale_color_manual(values = custom_pal) +
  labs(title = paste("selfunsuregodreal", "godvoxaloud", sep = " x "),
       # subtitle = "Excluding people who did not have a clear answer",
       x = "selfunsuregodreal",
       y = "godvoxaloud",
       color = "Researcher") +
  theme_bw() +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  guides(color = guide_legend(ncol = 6, byrow = TRUE, 
                              override.aes = list(alpha = 1)))

d %>%
  filter(!is.na(selfunsuregodreal), !is.na(godvoxaloud),
         selfunsuregodreal != "Other", godvoxaloud != "Other") %>%
  ggplot(aes(x = selfunsuregodreal, alpha = godvoxaloud, fill = researcher)) +
  facet_grid(urban_rural ~ charismatic_local ~ country) +
  geom_bar(position = "fill") +
  geom_text(data = selfunsuregodreal_count_by_quad,
            aes(x = selfunsuregodreal, y = 1, alpha = NULL, fill = NULL,
                label = paste0("(n=", n, ")")),
            size = 2, nudge_y = 0.05) +
  # scale_fill_brewer(guide = NULL, palette = "Dark2") +
  scale_fill_manual(guide = NULL, values = custom_pal) +
  scale_alpha_discrete(range = c(0.5, 1)) +
  scale_y_continuous(breaks = seq(0, 1, 0.25)) +
  labs(title = paste("selfunsuregodreal", "godvoxaloud", sep = " x "),
       subtitle = "Excluding people who did not have a clear answer",
       x = "selfunsuregodreal",
       y = "Proportion",
       alpha = "godvoxaloud",
       fill = "Researcher") +
  theme_bw() +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  guides(alpha = guide_legend(),
         fill = "none")

Ignoring unknown aesthetics: fill

d %>%
  filter(!is.na(selfunsuregodreal), !is.na(godvoxaloud),
         selfunsuregodreal != "Other", godvoxaloud != "Other") %>%
  ggplot(aes(x = selfunsuregodreal, alpha = godvoxaloud, fill = country)) +
  facet_grid(. ~ country) +
  geom_bar(position = "fill") +
  geom_text(data = selfunsuregodreal_count,
            aes(x = selfunsuregodreal, y = 1, alpha = NULL, fill = NULL,
                label = paste0("(n=", n, ")")),
            size = 2, nudge_y = 0.05) +
  scale_fill_brewer(guide = NULL, palette = "Dark2") +
  # scale_fill_manual(guide = NULL, values = custom_pal) +
  scale_alpha_discrete(range = c(0.5, 1)) +
  scale_y_continuous(breaks = seq(0, 1, 0.25)) +
  labs(title = paste("selfunsuregodreal", "godvoxaloud", sep = " x "),
       subtitle = "Excluding people who did not have a clear answer",
       x = "selfunsuregodreal",
       y = "Proportion",
       alpha = "godvoxaloud",
       fill = "Researcher") +
  theme_bw() +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  guides(alpha = guide_legend(),
         fill = "none")

Ignoring unknown aesthetics: fill

d4 <- d %>%
  filter(!is.na(selfunsuregodreal), !is.na(godvoxaloud),
         selfunsuregodreal != "Other", godvoxaloud != "Other",
         selfunsuregodreal != "Maybe") %>%
  mutate(selfunsuregodreal = factor(selfunsuregodreal,
                                    levels = c("No", "Yes")),
         selfunsuregodreal_num = as.numeric(selfunsuregodreal) - 1,
         godvoxaloud = factor(godvoxaloud,
                               levels = c("No", "Yes")),
         godvoxaloud_num = as.numeric(godvoxaloud) - 1) %>%
  select(country, researcher, urban_rural, charismatic_local, subject_name,
         starts_with("selfunsuregodreal"), starts_with("godvoxaloud"))
contrasts(d4$country) <- cbind("GH" = c(-1, 1, 0, 0, 0),
                               "TH" = c(-1, 0, 1, 0, 0),
                               "CH" = c(-1, 0, 0, 1, 0),
                               "VT" = c(-1, 0, 0, 0, 1))
contrasts(d4$selfunsuregodreal) <- cbind("Y" = c(-1, 1))
r4 <- glm(godvoxaloud_num ~ selfunsuregodreal + country,
                  family = "binomial", data = d4)
summary(r4)

Call:
glm(formula = godvoxaloud_num ~ selfunsuregodreal + country, 
    family = "binomial", data = d4)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.4270  -1.0266  -0.5155   0.9986   2.0415  

Coefficients:
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)        -0.47400    0.12983  -3.651 0.000261 ***
selfunsuregodrealY -0.06683    0.13151  -0.508 0.611315    
countryGH           0.97704    0.21634   4.516 6.29e-06 ***
countryTH           0.01075    0.26818   0.040 0.968028    
countryCH           0.38097    0.24572   1.550 0.121032    
countryVT           0.04141    0.26124   0.159 0.874046    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 426.65  on 313  degrees of freedom
Residual deviance: 381.81  on 308  degrees of freedom
AIC: 393.81

Number of Fisher Scoring iterations: 4

Spiritual experiences indexed by godviavisions

_Some people say that they have had a vision from God*—they have a picture, but it is like they see it with their eyes. Has anything like that happened to you?_

d %>%
  filter(!is.na(selfunsuregodreal), !is.na(godviavisions)) %>% #,
         # selfunsuregodreal != "Other", godviavisions != "Other") %>%
  ggplot(aes(x = selfunsuregodreal, y = godviavisions, color = researcher)) +
  facet_grid(urban_rural ~ charismatic_local ~ country) +
  geom_point(alpha = 0.3, position = position_jitter(0.2, 0.2)) +
  scale_color_manual(values = custom_pal) +
  labs(title = paste("selfunsuregodreal", "godviavisions", sep = " x "),
       # subtitle = "Excluding people who did not have a clear answer",
       x = "selfunsuregodreal",
       y = "godviavisions",
       color = "Researcher") +
  theme_bw() +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  guides(color = guide_legend(ncol = 6, byrow = TRUE, 
                              override.aes = list(alpha = 1)))

d %>%
  filter(!is.na(selfunsuregodreal), !is.na(godviavisions),
         selfunsuregodreal != "Other", godviavisions != "Other") %>%
  ggplot(aes(x = selfunsuregodreal, alpha = godviavisions, fill = researcher)) +
  facet_grid(urban_rural ~ charismatic_local ~ country) +
  geom_bar(position = "fill") +
  geom_text(data = selfunsuregodreal_count_by_quad,
            aes(x = selfunsuregodreal, y = 1, alpha = NULL, fill = NULL,
                label = paste0("(n=", n, ")")),
            size = 2, nudge_y = 0.05) +
  # scale_fill_brewer(guide = NULL, palette = "Dark2") +
  scale_fill_manual(guide = NULL, values = custom_pal) +
  scale_alpha_discrete(range = c(0.5, 1)) +
  scale_y_continuous(breaks = seq(0, 1, 0.25)) +
  labs(title = paste("selfunsuregodreal", "godviavisions", sep = " x "),
       subtitle = "Excluding people who did not have a clear answer",
       x = "selfunsuregodreal",
       y = "Proportion",
       alpha = "godviavisions",
       fill = "Researcher") +
  theme_bw() +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  guides(alpha = guide_legend(),
         fill = "none")

Ignoring unknown aesthetics: fill

d %>%
  filter(!is.na(selfunsuregodreal), !is.na(godviavisions),
         selfunsuregodreal != "Other", godviavisions != "Other") %>%
  ggplot(aes(x = selfunsuregodreal, alpha = godviavisions, fill = country)) +
  facet_grid(. ~ country) +
  geom_bar(position = "fill") +
  geom_text(data = selfunsuregodreal_count,
            aes(x = selfunsuregodreal, y = 1, alpha = NULL, fill = NULL,
                label = paste0("(n=", n, ")")),
            size = 2, nudge_y = 0.05) +
  scale_fill_brewer(guide = NULL, palette = "Dark2") +
  # scale_fill_manual(guide = NULL, values = custom_pal) +
  scale_alpha_discrete(range = c(0.5, 1)) +
  scale_y_continuous(breaks = seq(0, 1, 0.25)) +
  labs(title = paste("selfunsuregodreal", "godviavisions", sep = " x "),
       subtitle = "Excluding people who did not have a clear answer",
       x = "selfunsuregodreal",
       y = "Proportion",
       alpha = "godviavisions",
       fill = "Researcher") +
  theme_bw() +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  guides(alpha = guide_legend(),
         fill = "none")

Ignoring unknown aesthetics: fill

d5 <- d %>%
  filter(!is.na(selfunsuregodreal), !is.na(godviavisions),
         selfunsuregodreal != "Other", godviavisions != "Other",
         selfunsuregodreal != "Maybe") %>%
  mutate(selfunsuregodreal = factor(selfunsuregodreal,
                                    levels = c("No", "Yes")),
         selfunsuregodreal_num = as.numeric(selfunsuregodreal) - 1,
         godviavisions = factor(godviavisions,
                               levels = c("No", "Yes")),
         godviavisions_num = as.numeric(godviavisions) - 1) %>%
  select(country, researcher, urban_rural, charismatic_local, subject_name,
         starts_with("selfunsuregodreal"), starts_with("godviavisions"))
contrasts(d5$country) <- cbind("GH" = c(-1, 1, 0, 0, 0),
                               "TH" = c(-1, 0, 1, 0, 0),
                               "CH" = c(-1, 0, 0, 1, 0),
                               "VT" = c(-1, 0, 0, 0, 1))
contrasts(d5$selfunsuregodreal) <- cbind("Y" = c(-1, 1))
r5 <- glm(godviavisions_num ~ selfunsuregodreal + country,
                  family = "binomial", data = d5)
summary(r5)

Call:
glm(formula = godviavisions_num ~ selfunsuregodreal + country, 
    family = "binomial", data = d5)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.3530  -0.9346  -0.6884   1.0629   1.8214  

Coefficients:
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)        -0.61222    0.12652  -4.839 1.30e-06 ***
selfunsuregodrealY  0.06415    0.13016   0.493   0.6221    
countryGH           0.95182    0.21279   4.473 7.71e-06 ***
countryTH           0.33395    0.26614   1.255   0.2096    
countryCH          -0.58890    0.27649  -2.130   0.0332 *  
countryVT           0.07429    0.25358   0.293   0.7696    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 422.46  on 318  degrees of freedom
Residual deviance: 391.05  on 313  degrees of freedom
AIC: 403.05

Number of Fisher Scoring iterations: 4

Spiritual experiences indexed by godviabodyexperiences

Some people have particular experiences in your body that they associate with God* or spirit. Does that happen for you? [examples: warm hands, goosebumps, fire in the belly]

d %>%
  filter(!is.na(selfunsuregodreal), !is.na(godviabodyexperiences)) %>% #,
         # selfunsuregodreal != "Other", godviabodyexperiences != "Other") %>%
  ggplot(aes(x = selfunsuregodreal, y = godviabodyexperiences, color = researcher)) +
  facet_grid(urban_rural ~ charismatic_local ~ country) +
  geom_point(alpha = 0.3, position = position_jitter(0.2, 0.2)) +
  scale_color_manual(values = custom_pal) +
  labs(title = paste("selfunsuregodreal", "godviabodyexperiences", sep = " x "),
       # subtitle = "Excluding people who did not have a clear answer",
       x = "selfunsuregodreal",
       y = "godviabodyexperiences",
       color = "Researcher") +
  theme_bw() +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  guides(color = guide_legend(ncol = 6, byrow = TRUE, 
                              override.aes = list(alpha = 1)))

d %>%
  filter(!is.na(selfunsuregodreal), !is.na(godviabodyexperiences),
         selfunsuregodreal != "Other", godviabodyexperiences != "Other") %>%
  ggplot(aes(x = selfunsuregodreal, alpha = godviabodyexperiences, fill = researcher)) +
  facet_grid(urban_rural ~ charismatic_local ~ country) +
  geom_bar(position = "fill") +
  geom_text(data = selfunsuregodreal_count_by_quad,
            aes(x = selfunsuregodreal, y = 1, alpha = NULL, fill = NULL,
                label = paste0("(n=", n, ")")),
            size = 2, nudge_y = 0.05) +
  # scale_fill_brewer(guide = NULL, palette = "Dark2") +
  scale_fill_manual(guide = NULL, values = custom_pal) +
  scale_alpha_discrete(range = c(0.5, 1)) +
  scale_y_continuous(breaks = seq(0, 1, 0.25)) +
  labs(title = paste("selfunsuregodreal", "godviabodyexperiences", sep = " x "),
       subtitle = "Excluding people who did not have a clear answer",
       x = "selfunsuregodreal",
       y = "Proportion",
       alpha = "godviabodyexperiences",
       fill = "Researcher") +
  theme_bw() +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  guides(alpha = guide_legend(),
         fill = "none")

Ignoring unknown aesthetics: fill

d %>%
  filter(!is.na(selfunsuregodreal), !is.na(godviabodyexperiences),
         selfunsuregodreal != "Other", godviabodyexperiences != "Other") %>%
  ggplot(aes(x = selfunsuregodreal, alpha = godviabodyexperiences, fill = country)) +
  facet_grid(. ~ country) +
  geom_bar(position = "fill") +
  geom_text(data = selfunsuregodreal_count,
            aes(x = selfunsuregodreal, y = 1, alpha = NULL, fill = NULL,
                label = paste0("(n=", n, ")")),
            size = 2, nudge_y = 0.05) +
  scale_fill_brewer(guide = NULL, palette = "Dark2") +
  # scale_fill_manual(guide = NULL, values = custom_pal) +
  scale_alpha_discrete(range = c(0.5, 1)) +
  scale_y_continuous(breaks = seq(0, 1, 0.25)) +
  labs(title = paste("selfunsuregodreal", "godviabodyexperiences", sep = " x "),
       subtitle = "Excluding people who did not have a clear answer",
       x = "selfunsuregodreal",
       y = "Proportion",
       alpha = "godviabodyexperiences",
       fill = "Researcher") +
  theme_bw() +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  guides(alpha = guide_legend(),
         fill = "none")

Ignoring unknown aesthetics: fill

d6 <- d %>%
  filter(!is.na(selfunsuregodreal), !is.na(godviabodyexperiences),
         selfunsuregodreal != "Other", godviabodyexperiences != "Other",
         selfunsuregodreal != "Maybe") %>%
  mutate(selfunsuregodreal = factor(selfunsuregodreal,
                                    levels = c("No", "Yes")),
         selfunsuregodreal_num = as.numeric(selfunsuregodreal) - 1,
         godviabodyexperiences = factor(godviabodyexperiences,
                               levels = c("No", "Yes")),
         godviabodyexperiences_num = as.numeric(godviabodyexperiences) - 1) %>%
  select(country, researcher, urban_rural, charismatic_local, subject_name,
         starts_with("selfunsuregodreal"), starts_with("godviabodyexperiences"))
contrasts(d6$country) <- cbind("GH" = c(-1, 1, 0, 0, 0),
                               "TH" = c(-1, 0, 1, 0, 0),
                               "CH" = c(-1, 0, 0, 1, 0),
                               "VT" = c(-1, 0, 0, 0, 1))
contrasts(d6$selfunsuregodreal) <- cbind("Y" = c(-1, 1))
r6 <- glm(godviabodyexperiences_num ~ selfunsuregodreal + country,
                  family = "binomial", data = d6)
summary(r6)

Call:
glm(formula = godviabodyexperiences_num ~ selfunsuregodreal + 
    country, family = "binomial", data = d6)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.8606  -1.2321   0.6274   0.9309   1.2922  

Coefficients:
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)          0.7056     0.1277   5.523 3.33e-08 ***
selfunsuregodrealY   0.2421     0.1317   1.839   0.0659 .  
countryGH            0.5778     0.2363   2.446   0.0145 *  
countryTH            0.5884     0.3089   1.905   0.0568 .  
countryCH           -0.3357     0.2470  -1.359   0.1742    
countryVT           -0.1013     0.2519  -0.402   0.6876    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 411.83  on 319  degrees of freedom
Residual deviance: 390.07  on 314  degrees of freedom
AIC: 402.07

Number of Fisher Scoring iterations: 4

“Doubt” indxed by morequesmoreanswr

Do you think that the more spiritually mature you become, you will discover more questions or more answers?

morequesmoreanswr_count <- d %>%
  filter(!is.na(morequesmoreanswr), !is.na(godvoxaloud),
         morequesmoreanswr != "Other", godvoxaloud != "Other") %>%
  count(country, morequesmoreanswr) %>%
  data.frame()
morequesmoreanswr_count_by_quad <- d %>%
  filter(!is.na(morequesmoreanswr), !is.na(godvoxaloud),
         morequesmoreanswr != "Other", godvoxaloud != "Other") %>%
  count(country, urban_rural, charismatic_local, morequesmoreanswr) %>%
  data.frame()

Spiritual experiences indexed by godvoxaloud

Some people say that they have heard God* speak out loud to them. Has this ever happened to you?

d %>%
  filter(!is.na(morequesmoreanswr), !is.na(godvoxaloud)) %>% #,
         # morequesmoreanswr != "Other", godvoxaloud != "Other") %>%
  ggplot(aes(x = morequesmoreanswr, y = godvoxaloud, color = researcher)) +
  facet_grid(urban_rural ~ charismatic_local ~ country) +
  geom_point(alpha = 0.3, position = position_jitter(0.2, 0.2)) +
  scale_color_manual(values = custom_pal) +
  labs(title = paste("morequesmoreanswr", "godvoxaloud", sep = " x "),
       # subtitle = "Excluding people who did not have a clear answer",
       x = "morequesmoreanswr",
       y = "godvoxaloud",
       color = "Researcher") +
  theme_bw() +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  guides(color = guide_legend(ncol = 6, byrow = TRUE, 
                              override.aes = list(alpha = 1)))

d %>%
  filter(!is.na(morequesmoreanswr), !is.na(godvoxaloud),
         morequesmoreanswr != "Other", godvoxaloud != "Other") %>%
  ggplot(aes(x = morequesmoreanswr, alpha = godvoxaloud, fill = researcher)) +
  facet_grid(urban_rural ~ charismatic_local ~ country) +
  geom_bar(position = "fill") +
  geom_text(data = morequesmoreanswr_count_by_quad,
            aes(x = morequesmoreanswr, y = 1, alpha = NULL, fill = NULL,
                label = paste0("(n=", n, ")")),
            size = 2, nudge_y = 0.05) +
  # scale_fill_brewer(guide = NULL, palette = "Dark2") +
  scale_fill_manual(guide = NULL, values = custom_pal) +
  scale_alpha_discrete(range = c(0.5, 1)) +
  scale_y_continuous(breaks = seq(0, 1, 0.25)) +
  labs(title = paste("morequesmoreanswr", "godvoxaloud", sep = " x "),
       subtitle = "Excluding people who did not have a clear answer",
       x = "morequesmoreanswr",
       y = "Proportion",
       alpha = "godvoxaloud",
       fill = "Researcher") +
  theme_bw() +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  guides(alpha = guide_legend(),
         fill = "none")

Ignoring unknown aesthetics: fill

d %>%
  filter(!is.na(morequesmoreanswr), !is.na(godvoxaloud),
         morequesmoreanswr != "Other", godvoxaloud != "Other") %>%
  ggplot(aes(x = morequesmoreanswr, alpha = godvoxaloud, fill = country)) +
  facet_grid(. ~ country) +
  geom_bar(position = "fill") +
  geom_text(data = morequesmoreanswr_count,
            aes(x = morequesmoreanswr, y = 1, alpha = NULL, fill = NULL,
                label = paste0("(n=", n, ")")),
            size = 2, nudge_y = 0.05) +
  scale_fill_brewer(guide = NULL, palette = "Dark2") +
  # scale_fill_manual(guide = NULL, values = custom_pal) +
  scale_alpha_discrete(range = c(0.5, 1)) +
  scale_y_continuous(breaks = seq(0, 1, 0.25)) +
  labs(title = paste("morequesmoreanswr", "godvoxaloud", sep = " x "),
       subtitle = "Excluding people who did not have a clear answer",
       x = "morequesmoreanswr",
       y = "Proportion",
       alpha = "godvoxaloud",
       fill = "Researcher") +
  theme_bw() +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  guides(alpha = guide_legend(),
         fill = "none")

Ignoring unknown aesthetics: fill

d7 <- d %>%
  filter(!is.na(morequesmoreanswr), !is.na(godvoxaloud),
         morequesmoreanswr != "Other", godvoxaloud != "Other") %>%
  mutate(morequesmoreanswr = factor(morequesmoreanswr),
         morequesmoreanswr_num = as.numeric(morequesmoreanswr) - 1,
         godvoxaloud = factor(godvoxaloud,
                               levels = c("No", "Yes")),
         godvoxaloud_num = as.numeric(godvoxaloud) - 1) %>%
  select(country, researcher, urban_rural, charismatic_local, subject_name,
         starts_with("morequesmoreanswr"), starts_with("godvoxaloud"))
contrasts(d7$country) <- cbind("GH" = c(-1, 1, 0, 0, 0),
                               "TH" = c(-1, 0, 1, 0, 0),
                               "CH" = c(-1, 0, 0, 1, 0),
                               "VT" = c(-1, 0, 0, 0, 1))
contrasts(d7$morequesmoreanswr) <- cbind("Q" = c(1, -1))
r7 <- glm(godvoxaloud_num ~ morequesmoreanswr + country,
                  family = "binomial", data = d7)
summary(r7)

Call:
glm(formula = godvoxaloud_num ~ morequesmoreanswr + country, 
    family = "binomial", data = d7)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.4044  -0.9984  -0.4288   0.9900   2.2055  

Coefficients:
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)        -0.65574    0.16797  -3.904 9.46e-05 ***
morequesmoreanswrQ -0.03057    0.15819  -0.193    0.847    
countryGH           1.14447    0.23706   4.828 1.38e-06 ***
countryTH           0.18839    0.31250   0.603    0.547    
countryCH           0.28946    0.33757   0.857    0.391    
countryVT           0.03154    0.32558   0.097    0.923    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 325.90  on 241  degrees of freedom
Residual deviance: 280.75  on 236  degrees of freedom
AIC: 292.75

Number of Fisher Scoring iterations: 4

Spiritual experiences indexed by godviavisions

_Some people say that they have had a vision from God*—they have a picture, but it is like they see it with their eyes. Has anything like that happened to you?_

d %>%
  filter(!is.na(morequesmoreanswr), !is.na(godviavisions)) %>% #,
         # morequesmoreanswr != "Other", godviavisions != "Other") %>%
  ggplot(aes(x = morequesmoreanswr, y = godviavisions, color = researcher)) +
  facet_grid(urban_rural ~ charismatic_local ~ country) +
  geom_point(alpha = 0.3, position = position_jitter(0.2, 0.2)) +
  scale_color_manual(values = custom_pal) +
  labs(title = paste("morequesmoreanswr", "godviavisions", sep = " x "),
       # subtitle = "Excluding people who did not have a clear answer",
       x = "morequesmoreanswr",
       y = "godviavisions",
       color = "Researcher") +
  theme_bw() +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  guides(color = guide_legend(ncol = 6, byrow = TRUE, 
                              override.aes = list(alpha = 1)))

d %>%
  filter(!is.na(morequesmoreanswr), !is.na(godviavisions),
         morequesmoreanswr != "Other", godviavisions != "Other") %>%
  ggplot(aes(x = morequesmoreanswr, alpha = godviavisions, fill = researcher)) +
  facet_grid(urban_rural ~ charismatic_local ~ country) +
  geom_bar(position = "fill") +
  geom_text(data = morequesmoreanswr_count_by_quad,
            aes(x = morequesmoreanswr, y = 1, alpha = NULL, fill = NULL,
                label = paste0("(n=", n, ")")),
            size = 2, nudge_y = 0.05) +
  # scale_fill_brewer(guide = NULL, palette = "Dark2") +
  scale_fill_manual(guide = NULL, values = custom_pal) +
  scale_alpha_discrete(range = c(0.5, 1)) +
  scale_y_continuous(breaks = seq(0, 1, 0.25)) +
  labs(title = paste("morequesmoreanswr", "godviavisions", sep = " x "),
       subtitle = "Excluding people who did not have a clear answer",
       x = "morequesmoreanswr",
       y = "Proportion",
       alpha = "godviavisions",
       fill = "Researcher") +
  theme_bw() +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  guides(alpha = guide_legend(),
         fill = "none")

Ignoring unknown aesthetics: fill

d %>%
  filter(!is.na(morequesmoreanswr), !is.na(godviavisions),
         morequesmoreanswr != "Other", godviavisions != "Other") %>%
  ggplot(aes(x = morequesmoreanswr, alpha = godviavisions, fill = country)) +
  facet_grid(. ~ country) +
  geom_bar(position = "fill") +
  geom_text(data = morequesmoreanswr_count,
            aes(x = morequesmoreanswr, y = 1, alpha = NULL, fill = NULL,
                label = paste0("(n=", n, ")")),
            size = 2, nudge_y = 0.05) +
  scale_fill_brewer(guide = NULL, palette = "Dark2") +
  # scale_fill_manual(guide = NULL, values = custom_pal) +
  scale_alpha_discrete(range = c(0.5, 1)) +
  scale_y_continuous(breaks = seq(0, 1, 0.25)) +
  labs(title = paste("morequesmoreanswr", "godviavisions", sep = " x "),
       subtitle = "Excluding people who did not have a clear answer",
       x = "morequesmoreanswr",
       y = "Proportion",
       alpha = "godviavisions",
       fill = "Researcher") +
  theme_bw() +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  guides(alpha = guide_legend(),
         fill = "none")

Ignoring unknown aesthetics: fill

d8 <- d %>%
  filter(!is.na(morequesmoreanswr), !is.na(godviavisions),
         morequesmoreanswr != "Other", godviavisions != "Other") %>%
  mutate(morequesmoreanswr = factor(morequesmoreanswr),
         morequesmoreanswr_num = as.numeric(morequesmoreanswr) - 1,
         godviavisions = factor(godviavisions,
                               levels = c("No", "Yes")),
         godviavisions_num = as.numeric(godviavisions) - 1) %>%
  select(country, researcher, urban_rural, charismatic_local, subject_name,
         starts_with("morequesmoreanswr"), starts_with("godviavisions"))
contrasts(d8$country) <- cbind("GH" = c(-1, 1, 0, 0, 0),
                               "TH" = c(-1, 0, 1, 0, 0),
                               "CH" = c(-1, 0, 0, 1, 0),
                               "VT" = c(-1, 0, 0, 0, 1))
contrasts(d8$morequesmoreanswr) <- cbind("Q" = c(1, -1))
r8 <- glm(godviavisions_num ~ morequesmoreanswr + country,
                  family = "binomial", data = d8)
summary(r8)

Call:
glm(formula = godviavisions_num ~ morequesmoreanswr + country, 
    family = "binomial", data = d8)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.388  -1.039  -0.616   1.127   2.026  

Coefficients:
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)         -0.6647     0.1687  -3.939 8.18e-05 ***
morequesmoreanswrQ   0.1811     0.1525   1.188   0.2349    
countryGH            0.9661     0.2353   4.106 4.03e-05 ***
countryTH            0.5124     0.3171   1.616   0.1061    
countryCH           -1.0697     0.4511  -2.372   0.0177 *  
countryVT            0.3111     0.3162   0.984   0.3251    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 325.22  on 242  degrees of freedom
Residual deviance: 295.17  on 237  degrees of freedom
AIC: 307.17

Number of Fisher Scoring iterations: 4

Spiritual experiences indexed by godviabodyexperiences

Some people have particular experiences in your body that they associate with God* or spirit. Does that happen for you? [examples: warm hands, goosebumps, fire in the belly]

d %>%
  filter(!is.na(morequesmoreanswr), !is.na(godviabodyexperiences)) %>% #,
         # morequesmoreanswr != "Other", godviabodyexperiences != "Other") %>%
  ggplot(aes(x = morequesmoreanswr, y = godviabodyexperiences, color = researcher)) +
  facet_grid(urban_rural ~ charismatic_local ~ country) +
  geom_point(alpha = 0.3, position = position_jitter(0.2, 0.2)) +
  scale_color_manual(values = custom_pal) +
  labs(title = paste("morequesmoreanswr", "godviabodyexperiences", sep = " x "),
       # subtitle = "Excluding people who did not have a clear answer",
       x = "morequesmoreanswr",
       y = "godviabodyexperiences",
       color = "Researcher") +
  theme_bw() +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  guides(color = guide_legend(ncol = 6, byrow = TRUE, 
                              override.aes = list(alpha = 1)))

d %>%
  filter(!is.na(morequesmoreanswr), !is.na(godviabodyexperiences),
         morequesmoreanswr != "Other", godviabodyexperiences != "Other") %>%
  ggplot(aes(x = morequesmoreanswr, alpha = godviabodyexperiences, fill = researcher)) +
  facet_grid(urban_rural ~ charismatic_local ~ country) +
  geom_bar(position = "fill") +
  geom_text(data = morequesmoreanswr_count_by_quad,
            aes(x = morequesmoreanswr, y = 1, alpha = NULL, fill = NULL,
                label = paste0("(n=", n, ")")),
            size = 2, nudge_y = 0.05) +
  # scale_fill_brewer(guide = NULL, palette = "Dark2") +
  scale_fill_manual(guide = NULL, values = custom_pal) +
  scale_alpha_discrete(range = c(0.5, 1)) +
  scale_y_continuous(breaks = seq(0, 1, 0.25)) +
  labs(title = paste("morequesmoreanswr", "godviabodyexperiences", sep = " x "),
       subtitle = "Excluding people who did not have a clear answer",
       x = "morequesmoreanswr",
       y = "Proportion",
       alpha = "godviabodyexperiences",
       fill = "Researcher") +
  theme_bw() +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  guides(alpha = guide_legend(),
         fill = "none")

Ignoring unknown aesthetics: fill

d %>%
  filter(!is.na(morequesmoreanswr), !is.na(godviabodyexperiences),
         morequesmoreanswr != "Other", godviabodyexperiences != "Other") %>%
  ggplot(aes(x = morequesmoreanswr, alpha = godviabodyexperiences, fill = country)) +
  facet_grid(. ~ country) +
  geom_bar(position = "fill") +
  geom_text(data = morequesmoreanswr_count,
            aes(x = morequesmoreanswr, y = 1, alpha = NULL, fill = NULL,
                label = paste0("(n=", n, ")")),
            size = 2, nudge_y = 0.05) +
  scale_fill_brewer(guide = NULL, palette = "Dark2") +
  # scale_fill_manual(guide = NULL, values = custom_pal) +
  scale_alpha_discrete(range = c(0.5, 1)) +
  scale_y_continuous(breaks = seq(0, 1, 0.25)) +
  labs(title = paste("morequesmoreanswr", "godviabodyexperiences", sep = " x "),
       subtitle = "Excluding people who did not have a clear answer",
       x = "morequesmoreanswr",
       y = "Proportion",
       alpha = "godviabodyexperiences",
       fill = "Researcher") +
  theme_bw() +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  guides(alpha = guide_legend(),
         fill = "none")

Ignoring unknown aesthetics: fill

d9 <- d %>%
  filter(!is.na(morequesmoreanswr), !is.na(godviabodyexperiences),
         morequesmoreanswr != "Other", godviabodyexperiences != "Other") %>%
  mutate(morequesmoreanswr = factor(morequesmoreanswr),
         morequesmoreanswr_num = as.numeric(morequesmoreanswr) - 1,
         godviabodyexperiences = factor(godviabodyexperiences,
                               levels = c("No", "Yes")),
         godviabodyexperiences_num = as.numeric(godviabodyexperiences) - 1) %>%
  select(country, researcher, urban_rural, charismatic_local, subject_name,
         starts_with("morequesmoreanswr"), starts_with("godviabodyexperiences"))
contrasts(d9$country) <- cbind("GH" = c(-1, 1, 0, 0, 0),
                               "TH" = c(-1, 0, 1, 0, 0),
                               "CH" = c(-1, 0, 0, 1, 0),
                               "VT" = c(-1, 0, 0, 0, 1))
contrasts(d9$morequesmoreanswr) <- cbind("Q" = c(1, -1))
r9 <- glm(godviabodyexperiences_num ~ morequesmoreanswr + country,
                  family = "binomial", data = d9)
summary(r9)

Call:
glm(formula = godviabodyexperiences_num ~ morequesmoreanswr + 
    country, family = "binomial", data = d9)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9152  -1.2161   0.7310   0.9641   1.1674  

Coefficients:
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)         0.71079    0.15975   4.449 8.62e-06 ***
morequesmoreanswrQ  0.03348    0.15584   0.215   0.8299    
countryGH           0.43886    0.24850   1.766   0.0774 .  
countryTH           0.91566    0.38084   2.404   0.0162 *  
countryCH          -0.54853    0.33456  -1.640   0.1011    
countryVT          -0.15232    0.30627  -0.497   0.6190    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 310.96  on 244  degrees of freedom
Residual deviance: 294.81  on 239  degrees of freedom
AIC: 306.81

Number of Fisher Scoring iterations: 4

Absorption & spiritual experiences

Aggregate spiritual experiences score

d_spirit %>%
  filter(!is.na(abs_score), !is.na(spex_score)) %>%
  distinct(country, researcher, urban_rural, charismatic_local,
           subject_name, abs_score, spex_score) %>%
  ggplot(aes(x = abs_score, y = spex_score, color = researcher)) +
  facet_grid(urban_rural ~ charismatic_local ~ country) +
  geom_point(alpha = 0.3, position = position_jitter(0.2, 0)) +
  scale_color_manual(values = custom_pal) +
  labs(title = paste("abs_score", "spex_score", sep = " x "),
       # subtitle = "Excluding people who did not have a clear answer",
       x = "abs_score",
       y = "spex_score",
       color = "Researcher") +
  theme_bw() +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  guides(color = guide_legend(ncol = 6, byrow = TRUE, 
                              override.aes = list(alpha = 1)))

d_spirit %>%
  filter(!is.na(abs_score), !is.na(spex_score)) %>%
  distinct(country, researcher, urban_rural, charismatic_local,
           subject_name, abs_score, spex_score) %>%
  ggplot(aes(x = abs_score, y = spex_score, 
             color = country, fill = country, group = country)) +
  facet_grid(. ~ country) +
  geom_point(alpha = 0.3, position = position_jitter(0.2, 0)) +
  geom_smooth(aes(x = as.numeric(abs_score)), method = "lm") +
  # geom_smooth(aes(x = as.numeric(abs_score)), method = "loess") +
  # scale_color_manual(values = custom_pal) +
  # scale_fill_manual(values = custom_pal) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  labs(title = paste("abs_score", "spex_score", sep = " x "),
       # subtitle = "Excluding people who did not have a clear answer",
       x = "abs_score",
       y = "spex_score",
       color = "Researcher", fill = "Researcher") +
  theme_bw() +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  guides(color = guide_legend(ncol = 6, byrow = TRUE, 
                              override.aes = list(alpha = 1)),
         fill = guide_legend(ncol = 6, byrow = TRUE))

d_agg2 <- d_spirit %>%
  filter(!is.na(abs_score), !is.na(spex_score),
         abs_score != "Other", spex_score != "Other") %>%
  # distinct(country, researcher, urban_rural, charismatic_local,
  #          subject_name, abs_score, spex_score) %>%
  mutate(abs_score_num = as.numeric(abs_score) - 1) # %>%
  # select(country, researcher, urban_rural, charismatic_local, subject_name,
  #        starts_with("abs_score"), starts_with("spex_score"))
contrasts(d_agg2$country) <- cbind("GH" = c(-1, 1, 0, 0, 0),
                               "TH" = c(-1, 0, 1, 0, 0),
                               "CH" = c(-1, 0, 0, 1, 0),
                               "VT" = c(-1, 0, 0, 0, 1))
r_agg3 <- lmer(response ~ scale(abs_score, scale = F) + country 
               + (1 | subject_name) + (1 | question), 
               data = d_agg2)
summary(r_agg3)

Linear mixed model fit by REML ['lmerMod']
Formula: 
response ~ scale(abs_score, scale = F) + country + (1 | subject_name) +  
    (1 | question)
   Data: d_agg2

REML criterion at convergence: 8842.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.6193 -0.7848 -0.1847  0.8886  2.5512 

Random effects:
 Groups       Name        Variance Std.Dev.
 subject_name (Intercept) 0.02356  0.1535  
 question     (Intercept) 0.03631  0.1905  
 Residual                 0.17833  0.4223  
Number of obs: 7418, groups:  subject_name, 320; question, 24

Fixed effects:
                             Estimate Std. Error t value
(Intercept)                  0.446950   0.040202  11.118
scale(abs_score, scale = F)  0.005888   0.001420   4.146
countryGH                    0.028049   0.017994   1.559
countryTH                   -0.034780   0.021471  -1.620
countryCH                   -0.056734   0.021195  -2.677
countryVT                    0.143577   0.023114   6.212

Correlation of Fixed Effects:
            (Intr) s(_s=F cntrGH cntrTH cntrCH
scl(b_,s=F) -0.006                            
countryGH   -0.041 -0.051                     
countryTH    0.019 -0.012 -0.230              
countryCH    0.006  0.193 -0.219 -0.279       
countryVT    0.034 -0.242 -0.236 -0.298 -0.332

r_agg4 <- lmer(response ~ scale(abs_score, scale = F) * country 
               + (1 | subject_name) + (1 | question), 
               data = d_agg2)
summary(r_agg4)

Linear mixed model fit by REML ['lmerMod']
Formula: 
response ~ scale(abs_score, scale = F) * country + (1 | subject_name) +  
    (1 | question)
   Data: d_agg2

REML criterion at convergence: 8879.6

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.6200 -0.7842 -0.1839  0.8860  2.5509 

Random effects:
 Groups       Name        Variance Std.Dev.
 subject_name (Intercept) 0.02371  0.1540  
 question     (Intercept) 0.03630  0.1905  
 Residual                 0.17833  0.4223  
Number of obs: 7418, groups:  subject_name, 320; question, 24

Fixed effects:
                                        Estimate Std. Error t value
(Intercept)                            0.4455970  0.0403955  11.031
scale(abs_score, scale = F)            0.0064983  0.0015318   4.242
countryGH                              0.0301795  0.0185149   1.630
countryTH                             -0.0355340  0.0219127  -1.622
countryCH                             -0.0592259  0.0222275  -2.665
countryVT                              0.1443285  0.0263563   5.476
scale(abs_score, scale = F):countryGH -0.0015553  0.0024737  -0.629
scale(abs_score, scale = F):countryTH  0.0052519  0.0035136   1.495
scale(abs_score, scale = F):countryCH -0.0020379  0.0028943  -0.704
scale(abs_score, scale = F):countryVT -0.0004637  0.0035465  -0.131

Correlation of Fixed Effects:
            (Intr) sc(_,s=F) cntrGH cntrTH cntrCH cntrVT s(_,s=F):G
scl(b_,s=F) -0.026                                                 
countryGH   -0.060  0.020                                          
countryTH    0.001  0.012    -0.185                                
countryCH    0.006  0.206    -0.194 -0.258                         
countryVT    0.068 -0.308    -0.301 -0.333 -0.336                  
s(_,s=F):GH  0.006 -0.287    -0.103  0.012 -0.109  0.207           
s(_,s=F):TH  0.003  0.183     0.010 -0.067 -0.074  0.147 -0.258    
s(_,s=F):CH  0.060 -0.076    -0.111 -0.092  0.228  0.092 -0.129    
s(_,s=F):VT -0.087  0.196     0.205  0.176  0.089 -0.468 -0.265    
            s(_,s=F):T s(_,s=F):C
scl(b_,s=F)                      
countryGH                        
countryTH                        
countryCH                        
countryVT                        
s(_,s=F):GH                      
s(_,s=F):TH                      
s(_,s=F):CH -0.295               
s(_,s=F):VT -0.353     -0.299    

anova(r_agg3, r_agg4)

refitting model(s) with ML (instead of REML)

Data: d_agg2
Models:
r_agg3: response ~ scale(abs_score, scale = F) + country + (1 | subject_name) + 
r_agg3:     (1 | question)
r_agg4: response ~ scale(abs_score, scale = F) * country + (1 | subject_name) + 
r_agg4:     (1 | question)
       Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
r_agg3  9 8819.6 8881.8 -4400.8   8801.6                         
r_agg4 13 8825.1 8914.9 -4399.5   8799.1 2.5102      4     0.6428

d_r_agg3_predicted <- d_agg2 %>%
  mutate(response_pred = predict(r_agg3, d_agg2)) %>%
  group_by(subject_name) %>%
  mutate(spex_score_pred = sum(response_pred, na.rm = T)) %>%
  ungroup() %>%
  group_by(country, abs_score) %>%
  do(data.frame(rbind(smean.cl.boot(.$spex_score_pred))))
d_spirit %>%
  filter(!is.na(abs_score), !is.na(spex_score)) %>%
  distinct(country, researcher, urban_rural, charismatic_local,
           subject_name, abs_score, spex_score) %>%
  ggplot(aes(x = abs_score, y = spex_score, 
             color = country, fill = country, group = country)) +
  facet_grid(. ~ country) +
  geom_point(alpha = 0.3, position = position_jitter(0.2, 0)) +
  geom_line(data = d_r_agg3_predicted, aes(y = Mean)) +
  geom_ribbon(data = d_r_agg3_predicted, 
              aes(ymin = Lower, ymax = Upper, y = NULL), 
              alpha = 0.5, size = 0) +
  # scale_color_manual(values = custom_pal) +
  # scale_fill_manual(values = custom_pal) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  labs(title = paste("abs_score", "spex_score", sep = " x "),
       subtitle = "Lines are predictions from a mixed effects linear regression:\nlmer(response ~ abs_score + country + (1 | subject_name) + (1 | question)",
       x = "abs_score",
       y = "spex_score",
       color = "Site", fill = "Site") +
  theme_bw() +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  guides(color = "none",
         fill = "none")

Ignoring unknown aesthetics: y

Looking just at ‘apples’ (urban charismatics)

d_spirit %>%
  filter(!is.na(abs_score), !is.na(spex_score),
         quad == "urban charismatic") %>%
  distinct(country, researcher, charismatic_local, urban_rural, 
           subject_name, abs_score, spex_score) %>%
  ggplot(aes(x = abs_score, y = spex_score, 
             color = country, fill = country, group = country)) +
  facet_grid(charismatic_local ~ urban_rural ~ country) +
  geom_point(alpha = 0.3, position = position_jitter(0.2, 0)) +
  geom_smooth(method = "lm") +
  # scale_color_manual(values = custom_pal) +
  # scale_fill_manual(values = custom_pal) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  labs(title = paste("abs_score", "spex_score", sep = " x "),
       # subtitle = "Lines are predictions from a mixed effects linear regression:\nlmer(response ~ abs_score + country + (1 | subject_name) + (1 | question)",
       x = "abs_score",
       y = "spex_score",
       color = "Site", fill = "Site") +
  theme_bw() +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  guides(color = "none",
         fill = "none")

Looking at all charismatics

d_spirit %>%
  filter(!is.na(abs_score), !is.na(spex_score),
         charismatic_local == "charismatic") %>%
  distinct(country, researcher, charismatic_local, urban_rural, 
           subject_name, abs_score, spex_score) %>%
  ggplot(aes(x = abs_score, y = spex_score, 
             color = country, fill = country, group = country)) +
  facet_grid(charismatic_local ~ country) +
  geom_point(alpha = 0.3, position = position_jitter(0.2, 0)) +
  geom_smooth(method = "lm") +
  # scale_color_manual(values = custom_pal) +
  # scale_fill_manual(values = custom_pal) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  labs(title = paste("abs_score", "spex_score", sep = " x "),
       # subtitle = "Lines are predictions from a mixed effects linear regression:\nlmer(response ~ abs_score + country + (1 | subject_name) + (1 | question)",
       x = "abs_score",
       y = "spex_score",
       color = "Site", fill = "Site") +
  theme_bw() +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  guides(color = "none",
         fill = "none")

To do

Here are other things on our to-do list:

• Experiences by sense (Josh)
• Extremity of experiences (Josh)
• …